
    library(numDeriv)
    
    source("1.2.1 RDReg.R")

    CHR = function(data){
        # Input:
        #   data: a list containing
        #   (Y=Y, Z=Z, D = D, X=X, W=W, Delta=Delta, mono=mono)
        #   Y: observed survival time
        #   Z: instrument taking 0 or 1
        #   D: treatment taking 0 or 1
        #   X: Covariates used when mono = TRUE
        #   W: Covariates used when mono = FALSE
        #   Delta: censoring indicator, 1 = observed, 0 = censored
        
        MyAttach(data)
        
        if(is.vector(X)) p = 2 else p = ncol(X) + 1
        n = length(Y)
        
        pZ.m = glm(Z ~ X, family = binomial)
        pZ   = predict(pZ.m, type = "response")
        fZ = pZ; fZ[Z==0] = 1 - pZ[Z==0] 
        
        if(mono)    Va = cbind(rep(1,length(X)),X)        
        if(!mono)   Va = cbind(rep(1,length(W)),W)
        Vb = cbind(rep(1,length(X)),X)
        
        brm.rd = RDMLEst("RD", D, Z, Va, Vb)
        beta.hat  = brm.rd$point.est[1:p]; eta.hat = pZ.m$coefficients
        theta.hat = c(beta.hat, eta.hat)
        delta.D = tanh(Va %*% beta.hat)
        
        omega = (2*Z - 1) / (fZ * delta.D)
        gamma2 = sapply(1:n, function(i) sum(D * as.numeric(Y >= Y[i]) * omega))
        gamma1 = sapply(1:n, function(i) sum(as.numeric(Y >= Y[i]) * omega))
        
        e.psi = sum(Delta * D * omega * (gamma1 - gamma2)) /
            sum(Delta * (1-D) * omega * gamma2)
        
        # The point estimate for the causal hazard ratio
        psi.hat   = log(e.psi)
        
        if(is.na(psi.hat) | !is.finite(psi.hat)) return(c(psi.hat,NA))
        
        
        U = function(psi,theta){
            beta = theta[1:p];  eta = theta[(p+1):(2*p)]
            delta.D = tanh(Va %*% beta)
            V  = cbind(rep(1,length(X)),X)
            pZ = expit(V %*% eta)
            fZ = pZ; fZ[Z==0] = 1 - pZ[Z==0] 
            
            omega  = (2*Z - 1) / (fZ * delta.D)
            gamma2 = sapply(1:n, function(i) sum(D * as.numeric(Y >= Y[i]) * omega))
            gamma1 = sapply(1:n, function(i) sum(as.numeric(Y >= Y[i]) * omega))
            
            return(Delta * D * omega * (gamma1 - gamma2) * exp(-psi) -
                       (Delta * (1-D) * omega * gamma2))
        }
        mU = function(psi, theta) mean(U(psi,theta))
        mUpsi   = function(psi) mU(psi, theta.hat)
        mUtheta = function(theta) mU(psi.hat, theta)
        
        dmU.dpsi   = grad(mUpsi, psi.hat)
        dmU.dtheta = grad(mUtheta, theta.hat)
        IF.beta    = t(brm.rd$IF.alpha.beta[1:p,])
        
        V  = cbind(rep(1,length(X)),X)
        pZ = as.vector(expit(V %*% eta.hat))
        dpZ.deta = V * pZ * (1 - pZ)  # n by 2
        part4    = -t(V) %*% dpZ.deta  # 2 by 2
        pZ.score = V * (Z - pZ)
        IF.eta   =  - t((solve(part4)) %*% t(pZ.score))
        
        IF.theta   = cbind(IF.beta, IF.eta)
            
        Utilde     = U(psi.hat, theta.hat) + IF.theta %*% as.matrix(dmU.dtheta)
        
        IF.psi     = -(dmU.dpsi)^{-1} * Utilde
        
        # The standard deviation
        sd.psi    = sqrt(mean(IF.psi^2)/length(D))
            
        return(c(psi.hat, sd.psi))
        
    }